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Abstract 

We investigate the nature of a class of probability density functions, say G(a), with a the shape 
parameter, which generalizes the Gumbel distribution. These functions appear in a model of charge 
transport, when applied to a metal-insulator-metal structure, where the insulator is constituted by 
a monolayer of bacteriorhodopsin. Current shows a sharp increase above about 3 V, interpreted 
as the cross-over between direct and injection sequential-tunneling regimes. In particular, we show 
that, changing the bias value, the probability density function changes its look from bimodal to 
unimodal. Actually, the bimodal distributions can be resolved in at least a couple of G(a) functions 
with different values of the shape parameter. 

PACS numbers: 02.50.Ng Distribution theory and Monte Carlo studies 02.70.Rr General statistical methods 
05. 40. -a Fluctuation phenomena, random processes, noise, and Brownian motion 
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There is a widespread evidence that fluctuations of macroscopic observables exhibiting 



"extreme" events 



-yj can be described by means of the generalized Gumbel distribution 



G(a), with a the shape parameter: 

G(a) = 6 -^^exp{~a[e(a)(z + u(a)) + e- e ^ z+L '^}} (1) 
1 (a) 

where the function v(a) is defined in terms of the Gamma function T(a) and its derivatives 

v{a) 

r(a), ip{a) and 8 2 (a) indicating, respectively, the Gamma, digamma and trigamma function. 
In the definition given in (JTJ, G(a) is a normalized distribution function with zero mean and 
unitary variance. 

This distribution is quite intriguing, since it crosses the perimeter of the extreme value 
statistics, describing physical phenomena of very different nature . Among its main 

eatures, here we recall that the well-known Bramwell-Holdsworth-Pintor (BHP) distribution 

1 

7J, found relevant for the description of critical behaviors in magnetic, fluido dynamics 
and percolative systems [8[ is practically coincident with the distribution G(n/2) 
Furthermore, for a = 1 this distribution, also called scaled Gumbel distribution, is used to 
describe the fluctuations of the conditional galaxy density j9[. Finally, for a — > oo the G(a) 
recovers the Gaussian distribution. Therefore, this distribution has a unifying character since 
it describes, on respect the value of a, the statistics of events from critical to non-critical 
conditions. 

The shape parameter a can take all the positive real numbers jsj . The commonly accepted 
interpretation of integer valued G(a) is the straightforward generalization of the Gumbel 
distribution, i.e. it describes the distribution of the a-th largest value of a set of independent 
and identically distributed (iid) variables . Not so easy the interpretation of non-integer 
valued G(a): in this case distributions are related to long, but finite, range correlations which 
do not allow for a single variable to be dominant [5|. In other words, when a is not integer, 
G(a) describes a finite-size system with correlations at all the scale lengths. As a rough 
distinction, we can say that integer a values correspond to the establishment of ordering 
in the relevant variables, and non integer a values correspond to a loss of ordering, a more 
democratic situation in which variance is large and the weight of events far from the mean 
is large. 



In this letter we investigate the fluctuations around the steady state of the current flowing 
through a monolayer of bacteriorhodopsin (bR) [h]] , an integral membrane protein sensitive 
to the light, as function of the applied voltage. By making use of an atomic force microscope 
(AFM) technique, in a large range of applied bias (up to about 8 V), it was experimentally 



observed 



llj that the current exhibits a sharp transition between a near linear (Ohmic) 



regime and a superlinear one (roughly as V 8 increase of current at a threshold voltage of 
about 3 V). This behavior resembles a phase transition where the two different charge trans- 
port behaviors can be associated with a direct tunneling (DT) and an injection (or Fowler 
Nordheim) tunneling (FN) regime, respectively. This result was quantitatively reproduced 
within a microscopic model based on an impedance network protein analogous (INPA) 12 ]. 
Figured] reports the current voltage (I-V) characteristic, as obtained by the experiments 



ll| (see continuous curve) and the theoretical model (see the full circles) when the AFM tip 



just touches the protein monolayer at about 4.6 nm from the bottom metal-contact. The 
transport model uses a stochastic approach to select the tunneling mechanism (DT or FN) 
based on the probability reported in the inset of Fig. 1. 

By construction, the numerical approach allows for the simultaneous calculation of the 
current and its fluctuations around the steady state. Accordingly, with respect to the steady 
value, the calculated current evolution is found to exhibit spikes that resemble "extreme" 
events, whose number increases with the bias value, until it becomes difficult to establish 
what is e*t rem e and what nonnal fl 

The calculated variance of current fluctuations corresponding to the I-V characteristic of 
Fig. (JTJ, shows a a rather abrupt increase in concomitance with the cross-over region, at 
about 3 V. The giant increase, for about five orders in the magnitude of current variance, is 
associated with the opening of low resistance paths between contacts: they originate by the 
establishing of the FN regime which replaces the low voltage DT regime 13]. 

In the following we briefly survey the theoretical approach. The INPA model aims to 
predict the static and dynamical electrical responses of a protein in vitro, i.e. directly 
contacted or sandwiched with functionalized contacts to an external bias. In particular, the 
sing le protein is described by means of an irregular impedance network, with ideal contacts 



12J. The aminoacids constituting the tertiary structure of the protein are taken as the 



nodes of the network and the interaction between aminoacids as the links of the network. 
For a given applied voltage the network is solved within Kirchhoff rules by associating 
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FIG. 1: Experimental and calculated data for the I-V characteristics of Ref. In the inset, 
the transmission probability as given by DT (dashed line) and the interpolation of DT and FN 
(continuous line) for the typical parameters: m e the free electron mass, lij = 5.5 A, $=219 meV. 
The critical, Vc, and Ginzburg Vg, voltages are indicated. 



an impedance to each link. The transfer of electrons between a couple of nodes, selected 
according to an interaction radius, is ruled by two different probabilities: 



K = ex p 



for the DT regime and 



FN 
ij 



cxp 



■^J2m e ($-ie^.; 



2l iJ y/2m e \ $ 



eVj < $ , 



(2) 



h 



eV 



eVij > $ 



(3) 



for the FN regime. Here Vij is the local potential drop between the couple of i,j amino- 
acids and m e is the electron effective mass, here taken the same of the bare value. A smooth 
variation of the aminoacid resistivity is introduced to take into account the superlinear 
current response: 

Pmax eV < $ 

p(V) = I (4) 

where Pmax = 4 x 10 13 Vl A is the resistivity value which should be used to fit the I-V 
characteristic at the lowest voltages, p m %n = 4 x 10 5 Q A plays the role of an extremely low 
series resistance, limiting the current at the highest voltages, and $ = 219 meV is the value 
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FIG. 2: Histograms of current fluctuations, V=0.1-9V. 



of the energy barrier separating two nodes, here taken to be the same for all the couples of 



The electrical model well reproduces the measured I-V characteristic of bR. On the other 
hand, since the two tunneling regimes are stochastically chosen according to their probability, 
the model provides also the instantaneous current fluctuations. 

To estimate the probability density functions (PDFs) of current fluctuations, we follow 
this procedure. In a first step, we collect the histograms of ln(I) for different bias values 
in the range from 0.1 to 9 V. As reported in Fig. [2j in all the bias range the histograms 
strongly deviate from a symmetric Gaussian-like shape. Furthermore, at intermediate and 
high bias, i.e. near to the transition and beyond, they exhibit a nearly unimodal shape, 
while at bias lower than 1 V, the shape becomes bimodal. 

In a second step, we look for a fitting function of these histograms. For bias values larger 
than 1 V, the function: 



where A 0l Ai, a, m are the curve parameters, gives a good approximation of calculated data. 
To reduce the histogram to a standard PDF, the above parameters should be adjusted. 
In doing so, we guess that the final PDF is the generalized Gumbel distribution (CQ) with 
the same value of a; the standardization of histograms is performed by reshuffling the in- 
dependent variable as follows: x = az + b, and the new expression of the test function 



nodes [12 ] . 




(5) 
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H'(z) = A exp { -a {A x a){z + - + -) - a e -*"<*+£+?> 



becomes: 



which is proportional to G(a), when assuming: 

9. 



a a 



m 



(6) 



o = -± b= [ v a -->*. 

In particular, XH'(z) = G(a) and the frequency normalization is 

a a A x 
Y[a) A 

The parameters < x >, a, A are an a-generalization of the usual mean value (location 
parameter) and the scale parameter 14 1 . 



In the region of bias values l-j-5 V this procedure gives a single PDF, the "scaled Gumbel" 
G(l) [9j as shown in Fig. [3j 
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FIG. 3: Normalized distributions, V=l V (triangles), 2 V (empty squares), 5 V(diamonds). The 
continuous curve is G(l). 

In a third step, we analyze the bimodal histograms. Following the suggestions of 
the bimodal shape is interpreted as the effect of the superposition of different curves. In 
the present case we succeed in resolving two prevailing contributions: the i?(2)-histogram 
and a if(0.6)-histogram. By using the standardization technique shown above, both H(2) 
and if (0.6) can be traced back to the G(2) and G(0.6) distributions, respectively. The 
comparison between the Hs and Gs curves is shown in Figs. H] and El for bias V=0.1 V. 
Notice that the curve on G(2) could be also resolved in at least two different parts, thus 
signaling that the G(2) is a superposition of other PDFs. 
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FIG. 4: Histograms at V=0.1 V. Continuous line refers to the fitting function H(2), dashed line 
refers to the fitting function #(0.6), circles are the calculated data. 

In the last step we observe that far from the transition, the G(l) smooths down to a 
different unimodal distribution. Figure [6] reports the calculated data and the fit with two 
different PDFs. The G(0.7) gives the best fit. 
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FIG. 5: PDFs with V=0.1 V. Continuous line refers to the fitting function G(2), dashed line refers 
to the fitting function G(0.6), symbols are the calculated data when rescaled (see text). 



The most relevant result of the present investigation stems in solving the, apparently, 
substantial difference between the bimodal PDFs, evidenced before the transition, and the 
unimodal PDFs, emerging at the transition and beyond. Bimodal distributions have been 
previously investigated as "mixture" of universal unimodal distributions 15[, but to our 
knowledge, a continuous transformation from bimodal to unimodal behaviour has not been 



7 



-3 -2 -1 1 2 3 4 5 

(x-<x>)/a 



FIG. 6: PDFs of current fluctuations for V=9 V. Data have been compared with both the PDF 
(7(0.7) and G(0.6). To better resolve the differences among the distributions, the G(0.6) curve and 
the relative data are shifted of 0.1 on the vertical axis. 

observed. A further contribution to this topic is the observation that the PDF: 

G(a) = exp{-aw - ae w } (7) 

1 (a) 

is equivalent, in the sense of distribution, to the Gamma distribution of shape parameter a 
and life-time 1/A: 

f(t) = * A^V*. (8) 
This can be seen by using the change of variables: 

M- 



and requiring the normalization 



The Gamma distributions is well known in economics and granular materials [18] and 
is usually interpreted as the distribution of a sum of a (if a is integer) iid exponentially 
distributed variables characterized by the same life-time. In other words, the sum of different 
Gamma functions with the same life-time is still a Gamma function with a shape parameter 



equal to the sum of the single parameters 



17j |. Otherwise, for iid with different life times, 



the variab 



e sum (convolution) can be expressed by means of Gamma and Beta distribution 



functions 17j Therefore, we make the following conjecture: the G(2) PDF we found in the 



low bias region, takes into account two different kinds of current, one due to the maximal 



resistivity (DT regime) and the other one due to the minimal resistivity (FN regime) 13 ]. 
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The observed not-perfect convolution in a single G(2) could be due to the (very) different 
life times of the events. The G(l) PDF signals the existence of a single dominant current 
regime, due to the superposition of all the possible resistivities values (p(V), Pmax, Pmin) 
(FN). This perfect superposition is announced by the G(0.6) PDF which detaches by the 
G{2) at low bias. The shape parameter lower than 1 signals the not-complete superposition. 
In addition, above the transition, the minimal resistivity becomes more and more dominant 
and this drains the G(l) distribution producing a PDF with a shape parameter lower than 
1, G(0.7). 

The presence of G(l) in a wide bias range suggests a further interpretation of the phe- 
nomenon: the strengthening of the system configuration. It happens that the phase transi- 
tion between the two tunneling regimes is not abrupt but covers the bias value region 14- 5 
V. In this region, at the microscopic level, the links can choose among the resistivity values 
PMAXiJhfiin an d p{V) because both the DT and the FN regimes coexist. In different con- 



texts, [19] . the complete realization of the phase transition is associated with a temperature, 
called Ginzburg temperature, above which the ordering of the system is completed (in some 
cases with the formation of topological defects). By analogy, we can call the voltage value 
below which the transition is finalized (V=5 V), the Ginzburg voltage, Vq, and the voltage 
value at which the transition starts (Vc=l V), the critical voltage Vq (see Fig. [I]) [13| 
In conclusion, the main results presented in this letter are the following: 

1. The probability density functions of current fluctuations show both unimodal and 
bimodal shapes. 

2. The bimodal PDFs can be decomposed into at least two unimodal functions. 

3. All the PDFs can be drawn back to the parametric G(a) distribution function([T]). 

4. The shape parameter a is a function of the applied bias. 



As final remark, we observe that what found on current fluctuations is a model prediction, 
since no measurements have been performed so far. The proposed model finely describe 
the I-V characteristics, and the interpretation of the PDF of fluctuations is in line with 
the theoretical mechanism of current transport leaving the experimental test as a future 
challenge. 
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